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Abstract 

We consider finite element discretizations of the Biot’s consolidation model 
in poroelasticity with MINI and stabilized Pl-Pl elements. We analyze the 
convergence of the fully discrete model based on spatial discretization with 
these types of finite elements and implicit Euler method in time. We also 
address the issue related to the presence of non-physical oscillations in the 
pressure approximation for low permeabilities and/or small time steps. We 
show that even in ID a Stokes-stable finite element pair fails to provide a 
monotone discretization for the pressure in such regimes. We then introduce 
a stabilization term which removes the oscillations. We present numerical 
results confirming the monotone behavior of the stabilized schemes. 

Keywords: Stable finite elements, monotone discretizations, poroelasticity. 


1. Introduction 

The theory of poroelasticity models the interaction between the deforma¬ 
tion and the fluid flow in a fluid-saturated porous medium. Such coupling 
was already modelled in the early one-dimensional work of Terzaghi, see [T] , 
whereas the general three-dimensional mathematical model was established 
by Maurice Biot in several pioneering publications (see [2j and 0). 
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We assume here that the porous medium is linearly elastic, homoge¬ 
neous, isotropic and saturated by an incompressible Newtonian fluid. Un¬ 
der these assumptions, the quasi-static Biot’s model can be written as a 
time-dependent system of partial differential equations in the variables of 
displacements of the solid, u, and pressure of the fluid, p, 


— div a + Vp = /, a = 2pe(u) + A cliv(w)/ 

— div it + div K Vp = g, 


( 1 ) 

( 2 ) 


where a and £ are the effective stress and strain tensors, A and /j are the Lame 
coefficients, K is the hydraulic conductivity tensor, the right-hand term / is 
the density of applied body forces and the source term g represents a forced 
fluid extraction or injection process. The time derivative of the displacement 
vector is denoted by it. Results on the existence and uniqueness of the 
solution for these models have been investigated by Showalter in [3] and by 
Zenisek in [5], and the well-posedness for nonlinear poroelastic models is 
considered, for example, in [B]. 

Biot’s models are still used today in a great variety of fields, ranging 
from geomechanics and petroleum engineering, where these models have 
been applied ever since their discovery, to biomechanics or even food pro¬ 
cessing more recently. Some examples of applications in geosciences include 
petroleum production, solid waste disposal, carbon sequestration, soil consol¬ 
idation, glaciers dynamics, subsidence, liquefaction and hydraulic fracturing, 
for instance. In biomechanics the poroelastic theory can be used to describe 
tumor-induced stresses in the brain (see H). which can cause deformation 
of the surrounding tissue, and bone deformation under a mechanical load 
(see 0), for example. More recently, a promising and innovative application 
studies the food processes as a multiphase deformable porous media, in order 
to improve the quality and safety of the food, see (9j. 

Although some analytical solutions have been derived for some linear 
poroelasticity problems, see [TO] , and even some of them are obtained arti¬ 
ficially as in HU, numerical simulations seem to be the only way to obtain 
quantitative results for real applications. The numerical solution of these 
problems is usually based on finite element methods, see for example the 
monograph of Lewis and Schrefler in ra and the papers in [T3| HU .15] HE] . 
Finite difference methods have been also applied to solve this problem, see 
for example the convergence analysis in na and the extension to the discon¬ 
tinuous coefficients case in [IB. [Tflj. 
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It is well-known that approximations by standard finite difference and fi¬ 
nite element methods of the poroelasticity equations often exhibit strong non¬ 
physical oscillations in the fluid pressure, see for instance [III ED Eg ESI Eg. 
For example, this is the case when linear finite elements are used to approx¬ 
imate both displacement and pressure unknowns, or when a central finite 
difference scheme on collocated grids is considered. To eliminate such insta¬ 
bilities, approximation spaces for the vector and scalar fields, satisfying an 
appropriate inf-sup condition (see E3) are commonly used. Such discretiza¬ 
tions have been theoretically investigated by Murad et al. in [26J E3 l28j . 
As we show later, however, an inf-sup stable pair of spaces does not neces¬ 
sarily provide oscillation-free solutions. On the other hand, the oscillations 
disappear on very fine grids, but evidently, this is not always practical. 

Our work here is on investigating mechanisms for avoiding the nonphys¬ 
ical oscillations in the discrete solution, for example, by adding stabilization 
terms to the Galerkin formulation, while still maintaining the accuracy of 
approximations. Such strategy has been applied in [29j to provide a stable 
scheme by using linear finite element approximations for both unknowns. 
This was accomplished by adding an artificial term, namely, the time deriva¬ 
tive of a diffusion operator multiplied by a stabilization parameter, to the flow 
equation. The stabilization parameter, which depends on the elastic proper¬ 
ties of the solid and on the characteristic mesh size, was given a priori, and 
its optimality was shown in the one-dimensional case. This scheme provided 
solutions without oscillations independently of the chosen discretization pa¬ 
rameters. 

In this work, we present convergence analysis of fully discrete implicit 
schemes for the numerical solution of Biot’s consolidation model. We derive 
appropriate stabilization terms for both MINI element and Pl-Pl discretiza¬ 
tions, and numerically show that such choices of stabilization parameters and 
operators remove the non-physical oscillations in the approximations of the 
pressure. In this regard, our work fills in a gap in the literature, since to our 
knowledge the results presented here are the first theoretical results for fully 
discrete schemes involving stabilized spatial discretizations aimed to improve 
the monotonicity properties of the finite element schemes. 

The rest of the paper is organized as follows. In Section [2j we provide one 
dimensional example elements illustrating the undesirable oscillatory pres¬ 
sure behavior. We show both numerically and theoretically, that adding 
appropriate stabilization terms provide monotone discrete schemes and we 
calculate the exact values of the optimal stabilization parameters for both 
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MINI and Pl-Pl schemes. In Section |3] we show several abstract results on 
stabilized discretizations which we use in Section [4]to analyze the convergence 
of the fully discrete model. The abstract results in Section [3] apply to more 
general saddle-point problems with stabilization terms. In this section, we 
have also computed the exact Schur complement corresponding to the bubble 
functions in the MINI element. Next, in Section[4]we use the abstract results 
and show first order convergence in time and space for the fully discrete Biot’s 
consolidation model. The section [5] is devoted to the numerical study of the 
convergence and monotonicity properties of the resulting discretizations. We 
use several benchmark tests in poromechanics and show that appropriate 
choice of stabilization parameters result in approximations which respect the 
underlying physical behavior and are oscillation-free. Conclusions are drawn 
in Section [6l 


2. Pressure oscillatory behaviour: one dimensional example 


We consider an example modeling a column of height H of a porous 
medium saturated by an incompressible fluid, bounded by impermeable and 
rigid lateral walls and bottom, and supporting a load a 0 on the top which is 
free to drain. We have the following PDEs describing this model: 


d f du\ dp 

dx \ dx) + dx 
d fdu\ d / dp\ 

dt \dx) dx \ dx) 


(x,t) e (0 ,H) x (0, T\, 

0, 


( 3 ) 


with boundary and initial conditions 

du 

= p(o,*) = o, te(o,T], 

u(H,t)= 0, K^(H,t) = 0, fe(0,T], 

Ou 

— (a:,0) =0,i6 [0 , H], 


where E is the Young’s modulus and K is the hydraulic conductivity. It 
can be easily seen that problem (J3| is decoupled, giving rise to the following 
heat-type equation for the pressure 


d ( 1 \ d f dp\ 
dt\E P ) dx\ K dx) 


0 . 


( 4 ) 
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In order to discretize problem (|3|, we consider a non-uniform partition of 
spatial domain 17 = (0, H), 

0 = xo < x\ < ... < x n -\ < x n = H. 

In this way, the domain 17 is given by the disjoint union of elements Tj = 
[xi,Xi- (_i], 0 < i < n — 1, of size hi = Xi + \ — Xi. We assume that the Young 
modulus E(x) and the hydraulic conductivity K(x) are constants Ei and Ki 
on each element T{. Next, we are going to analyze two discretizations by two 
different pairs of finite elements with a backward Euler method in time. 


2.1. Discretization with linear finite elements 

First, we discretize using linear finite elements for both displacement and 
pressure. In this case, the following linear system of equations has to be 
solved on each time step 


r Ai 

Gi 


r up i 


0 

0 ' 


r urn-1 I 


k 

t A p _ 


L 

pm 

— 

i 

0 


L 

pm— 1 

+ 


(5) 


where m > 1, and r is the time discretization parameter. It is clear that the 
pressure at time level m must satisfy the following equation 


(Ci + rA p )P m = CiP m - 1 - GjA-'ifr - /r 1 ), (6) 

where Ci = —GjA^ l Gi is a tridiagonal matrix such that for an interior 
node Xi it is given by 


(c [ p m )i = i (Yi/r i + 


hi -1 lfih \ pm I pm 

7? ' W I ^ ' TP 1 

X'i—l WiJ XLi 


(7) 


Notice that the scheme associated with the above equation should be an ap¬ 
propriate discretization for problem Q. Depending on the relation between 
the space and time discretization parameters, the off-diagonal elements of 
matrix Ci + rA p could be positive and therefore the cause of possible non¬ 
physical oscillations in the approximation of the pressure. To avoid these 
instabilities, the following restriction holds, 


Ci 


max — ^ 
0<i<n—l AK.;E, 


< r. 


( 8 ) 


For example, in the case of an uniform-grid of size h and constant values 
of the parameters E and K in the whole domain, such restriction becomes 
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Figure 2.1: Numerical solution for the pressure field obtained with finite 
elements Pl-Pl and corresponding exact solution for (a) h = 1/32 and (b) 
h = 1/500. 


h 2 < AEKt. To confirm these unstable behavior, we solve system (J3]) in 
the computational domain (0,1) by using linear finite elements considering 
K Et = 10 -6 . In this case, it is necessary a mesh of at least 500 nodes to 
fulfill the restriction. In Figure [2T] we show the corresponding approximation 
of the pressure at the first time step, for two different values of h, that is, 
(a) h = 1/32 and (b) h = 1/500. Besides, we have plotted the analytical 
solution of the problem (see m- We can observe that strong non-physical 
oscillations appear for this type of finite element approximations, when the 
space discretization parameter is not small enough. It is clear that this is 
due to a lack of monotonicity of the scheme. At a first glance, it appears 
that these oscillations might be related to the locking effect and/or the fact 
that the pair of finite element does not satisfy an inf-sup condition. However, 
since our test is an one-dimensional problem, elastic locking can not appear, 
and therefore, in general, this can not be the only cause of this oscillatory 
behavior. 

2.2. Discretization with Taylor-Hood elements 

We consider the Taylor-Hood finite element method proposed in |30j ap¬ 
proximating the displacement by continuous piecewise quadratic functions 
and the pressure by continuous piecewise linear functions. It is well-known 
that this pair of finite elements provides a stable discretization for the Stokes 
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equation and satisfies inf-sup condition. Following similar computations as 
for the Pl-Pl case, and we obtain the following linear system of equations 
on each time step 


A b 0 Gb 


~u?~ 


i 

O 

O 

o 


r ur 1 1 


-pm 

Jb 

0 Ai Gi 


ur 

= 

0 0 0 


ur 1 

+ 

pm 

Jl 

G T b Gf tA p _ 


pm 


_G T b Gf 0_ 


pm— 1 


0 


where Ai, Gi correspond again to the linear basis functions whereas A b , Gb 
are associated with the bubble basis functions. In this case, the pressure at 
time level m satisfies the equation 


(c,+a+r^)p’" = (Q+a)p m - 1 -Gf^r 1 (/r-/r 1 )-G^ 1 (/r-/r 1 ), 

( 10 ) 

where Ci is as in ([7]) and Cb = —G^A b is given by 


(C b P m )i 


1 

12 


^*-1 pm | ( hi- 1 hi \ 

Ej « 



Note that the off-diagonal entries of matrix C b are non-positive, but again 
depending on the values of the parameters, the whole matrix Ci + C b + rA p 
can still have positive off-diagonal terms. To avoid this, on each element the 
restriction 


max 


0<i<n-l 6KiEi 


< r. 


( 11 ) 


must be fulfilled. 

In summary, the use of quadratic finite elements for displacement does 
contributes towards the reduction of the non-physical oscillations, but is still 
not enough to eliminate them. 

To illustrate this behavior, we consider again system 0 on an uniform 
grid of size h and constant coefficients E and K. In this particular case, 
the restriction ( [II] ) is simplified to h 2 < 6 EKr, and when EKr = 10 6 it 
is deduced that 409 nodes are needed to ensure a non-oscillatory behavior. 
In Figure ( 2.2[ ) we show the corresponding approximation of the pressure 
at the first time step, for two different values of h, that is, h = 1/32 and 
h = 1/409. Notice again that in the first case the pressure is not monotone 
(oscillations show up), which shows that the inf-sup condition is not enough 
for the monotonicity of the discretization. 
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Figure 2.2: Numerical solution for the pressure field obtained with finite 
elements P2-P1 and corresponding exact solution for (a) h = 1/32 and (b) 
h = 1/409. 


2.3. Monotone discretizations using perturbations 

To avoid the restrictions (|8]) for Pl-Pl and 0 for P2-P1 which result 
in the requirement for using very small mesh size, we are going to introduce 
a perturbation which will lead to monotone (and accurate) discretization 
independently of the chosen parameters. 

One way to achieve this is to add stabilization terms so that the dis¬ 
cretizations (J6]) and (10) correspond to the standard monotone linear finite 
element discretization of the parabolic (heat) equation Q. We define the 
following tridiagonal matrix 


(A e P m )i = £ 


hi —1 

Ei -1 


pm , 
- 0-1 ^ 


hi -1 , hi \ _ hi_ pm 

rp ' tp ) * tp O+i 

O-i O / o 


( 12 ) 


where e = 1/4 for the linear finite element pair and £ = 1/6 for the Taylor- 
Hood method. Then, it is clear that the perturbation of scheme ([6]) 


(Q + A e + rA p )P m = (Cj + A e )P m ~ 1 - Gj A^if? ~ fi ), 


(13) 


or the perturbation of (10) 


(' C l +C b +A e +rA p )P m = (Q+a+A^P^-Gi Ay (fJ n —fl n )—G{ ) 

(14) 
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Figure 2.3: Numerical solution for the pressure field obtained with the sta¬ 
bilized finite elements (a) Pl-Pl and (b) P2-P1 and corresponding exact 
solution. 


gives the standard discretization of Q by linear finite element method with 
mass-lumping. We also note that this perturbation corresponds to adding 
the following term to the second equation in ([3]) 


n— 1 





/vryvpfj 


V<7h dx. 


(15) 


Finally, in Figure [273] we show the approximation for the pressure obtained 
using the stabilized scheme for both the linear finite element pair and the 
Taylor-Hood method with h = 1/32 and we obtain monotone approximation 
for the pressure. 


3. Stability of discretizations and perturbations of Biot’s model 

In this section we provide results on the stability of discretizations of 
saddle point problems that can be viewed as perturbations of the Stokes 
equations. By stability, here, we mean bounds on the inverse of the discrete 
operator (for a fixed time step). We prove inf-sup condition for different 
discretizations for the poroelasticity problem, more precisely for MINI ele¬ 
ment and stabilized Pl-Pl schemes. Such results are well-known for Stokes 
equations (see, e.g. mmm)- 
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We hope that the results given below in Section 3.1 will be useful in other 
situations. We note that the generality of the abstract results allows us to use 
an unweighted L 2 norm for the pressure (not only an energy norm), which 
gives new estimates in the analysis of the fully discretized time dependent 
Biot’s model. 


3.1. Stability of a class of saddle point problems with perturbation 
In this section, we consider operators of the form 

A c =(i B ' c y.VxQ^V'x Q\ (16) 

where V and Q are Hilbert spaces and V and O' are their dual spaces. Here, 
(-, •) is the standard duality pairing and B' : Q H > V' is the adjoint of B. We 
make the following assumptions on A and C. 

(Al) The operator A : V H» V' is bounded, selfadjoint and positive definite. 
Thus, A provides a scalar product (•, -)a = (A-, •) and a norm on V 
denoted by || • ||^. The Hilbert Space V is then equipped with this 
inner product and norm, and we have that 

H 2 a :=(Av,v), \\f\\v> ■— (fiA~ l f), for all v € V, / G V 

II^IIvm-v' = ||"4 _1 ||v , i-^u = 1- 

(A2) The operator B : V i-)- Q' is bounded. 

(A3) Similarly to A, the operator C : Q ^ Q' is bounded, selfadjoint and 
positive (semi)definite. Thus on Q we have a norm (or a semi-norm) 
denoted by || ■ ||c 

We introduce a norm on V x Q: 

\\\(u,p)\i 2 =\\u\\ 2 A +\\p\\ 2 c +\\p\\ 2 . (17) 

We note that if C is only semidefinite, then || ■ ||c is only a seminorm on Q. 
Here || • || denotes the norm on Q and ||| • ||| is the norm on V x Q in which 
we will prove stability estimates for the operator Ac- 

Clearly, Ac can be viewed as a perturbation of Ao, he. the operator 
with (7 = 0. For detailed discussion on perturbations of such saddle point 
problems, we refer the reader to the recent monograph by Bofli, Brezzi and 
Fortin [33] , 
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We now state and prove a necessary and sufficient condition for Ac to 
be isomorphism under the assumptions (Al)-(A3). More general results also 
hold (with A only invertible on a subspace, etc), but to prove them would 
require more elaborate arguments and such generality is beyond the scope of 
our considerations here. We have the following theorem. 


Theorem 1. Assume that (Al)-(A3) hold. Then Ac defined in (16) is 
an isomorphism if and only if the operator B satisfies the following inf-sup 
condition: For any q G Q we have 


sup 

vev 


{Bv, q) 

IMU 


> 7 b|M| - |M|c 


(18) 


Proof. We first assume that (18) holds and we introduce the bilinear form 


(Ac{u,p)\ (v,q)) = ( Au,v) + (Bv,p) + (Bu,q) - (Cp, q) 


ft is easy to verify that the operator Ac is bounded in ||| • ||| since both A 
and B are continuous. From the inf-sup condition (18), for any p, there exist 
w E V, such that (Bw,p) > (7s||p|| — ||p||c)||HU- Since this inequality 
does not change when we multiply w by a positive scalar, without loss of 
generality, we may assume that ||u;||a = ||p||. We then have, 


(Bw,p) > (tbIMI - IHIc)IHI- 


For a given pair (u,p) E V x Q and with w defined as above, we choose 
v = u + Ow, and, q = —p, with some 6 > 0 to be determined later. Using 
the inf-sup condition, the fact that ||u;||^ = ||p|| and applying some obvious 
inequalities, such as, ab > —^a? — f& 2 , we have 

(Ac(u,p)] (v,q)) = (Au, u + Ow) + (B(u + 9w),p) - (Bu,p) + (Cp,p) 

— \\u\\ 2 a + 6{Au,w) + 6(Bw,p) + \\p\\c 

> ^IMIa - j\\p\\ 2 + »7 S ||p|| 2 - »||p|| c ||p|| + IIpIIc 

> i|l-ll 2 4 + («7 b - y) llPir - e (^iipIIJ + ^||p|| 2 ) + HpIIJ. 

Since the inequality above holds for any 6 > 0, we choose 6 — to obtain 
that 

(A c (u,p)](v,q)) > ^|HlU^N 2 + ^lbl| 2 c>7lll(^p)f 
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where 7 = 1 min{2, 7 ^}. On the other hand, the triangle inequality implies 
that 

I(u, q )I = I(u + 9 w,p )I <71 III (u,p) I, 

with 71 depending only on 7 #. Hence, 


(A c (u,p)] {■v,q )) 

IIIK?) II 


>7||(«,P)I, 


7 


1_ 

7i 


which shows that Ac is an isomorphism. 

To prove the other direction, that the invertibility of Ac implies condi- 

- - a — a 1/ Vinoo A „ I 


tion (18), for any q £ Q, we define v q = —A 1 B'q E V. Since Tic * Uq ' — 
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0 

Bv q — Cq 


the invertibility of Tic implies that 


Ml < IK,g)|| < Me II \\ B v q -Cq\\ Q , < \\A C II (\\Bv q \\ Q/ + \\Cq\\ Q/ ). 

Since C is symmetric and positive (semi)-dehnite, we have {Cq, s) < y/ {Cq, q) y/ {Cs, s ). 
Hence, 

IIMIIq' = sup < y/\\C\\{Cq,q). 

seQ ||s|| 

To estimate ||I?u 9 ||q/ we observe that ||Ht> g ||Q/ = sup sg g and we also 

have for all s G Q, 

\{Bv q ,s)\ = \{B's, A~ 1 B'q)\ < \{B's, A^B'q)\ 


B's 


v> 


^ llD/|l (/> A~ l B’q) ||D/|| {Aw,A~'B'q) 

f&V' \\f\\v wev ||Hwi||v" 

. lM ( Bw,q) ...... ( Bw,q) 

< ||H || ||H _1 || sup ,, 7 = ll-B || sup . 


uev F i 


?ev Fa 


The inf-sup condition (18) easily follows by combining the last two estimates. 

□ 


We have the following immediate corollaries. 

Corollary 1. Suppose that (Al)-(A3) hold. If Aq is an isomorphism, then 
Tic is an isomorphism for all continuous and positive (semi-) definite C. 
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Proof. From the fact that Aq is isomorphism it follows that (18) holds with 


C — 0, and hence, also with any symmetric positive, (semi-)dehnite and 
bounded C. This in turn (by Theorem [I]) implies that Ac is an isomorphism. 

□ 

The next corollary allows us to add consistent perturbations to already 
stable discretizations in order to improve the monotonicity properties of the 
underlying discretizations. 

Corollary 2 . Suppose that Ac is an isomorphism, that (Al)-(A3) hold, and 
that D is spectrally equivalent to C, namely a 0 |kl|c < \\q\\d < a ilkllc f or 
some positive constants «o and ol Then Ad is an isomorphism. 


Proof. For all q G Q, we have 
(Bv, q) 


d 


■ sup 

v£V 


> min{l, a 0 } 


c + sup 
vev 


( Bv , q) 


> min{l, «o}7_b 11 ^ 11 , 


\ v \\a V v&v ||U|U 

which shows (18) for Ad- Applying Theorem [I] gives the desired result. Q 


3.2. Application to discretizations of Biot’s model 

After a time discretization (backward Euler scheme in time) of the Biot’s 
model, the following system of differential equations is solved on every time 
step on a domain C M. d : 

— diva + Vp =/, a = 2pe(u) T Adiv(-u)/ (19) 

— div u + t div KVp = g. (20) 

A typical set of boundary conditions is 

u = 0, and (A'Vp • n ) = 0, on r c , 
a ■ n — j3 , and p = 0 on r f . 

To introduce the spatial discretization of the Biot’s model, we consider 
Enite dimensional spaces 14 C [f/p c (f2)] d and Qh C Hf t (Q) where Hf c (Q) 
and Hf ( fl) are the standard Sobolev spaces with functions whose traces 
vanish on r c and F t respectively. 

We have the following discrete formulation (on each time step) corre¬ 


sponding to ([19])—([20]) . Find (u,p) e 14 x Q h such that 

for all v G 14 , 
for all q G Qh- 


a(u,v) - (div v,p) = (f,v), 

— (div u, q) - ra p (p , q) = (g, q ), 


( 21 ) 

( 22 ) 
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The bilinear form a(-, •) is as follows: 


a(u,v ) = 2 n e{u) : e{v) + A / div u div w, a p (p,q) = / A'Vp ■ X7q. 

Jo, Jci J cl 

The corresponding operators A : 14 ^ 14', B : Qh ^ V^, and the norm on 
Q , || • || 5 are defined as follows: 

(Am,v) :=a(u,v), {Bu, q) := -(divu,g), (A p p, q) :=a p (p,q), 

\\q\\ 2 -.= T(A p q,q) + \\q\\l 2m . 

Since C may take different form for different discretizations, we do not specify 
its definition here. 


3.2.1. Discretization with MINI element 

We consider a discretization with MINI element, introduced in j3T] where 
the finite element spaces that we use are as follows: 

V h x Q h , where V h = V t ® V b , 

where VJ is the space of piece-wise (with respect to a triangulation 7/J lin¬ 
ear continuous vector valued functions on and V b is the space of bubble 
functions, defined as 

V b = span{<p b>T ei,..., ^p b ,T^d}TeT h -> Vb,T = ■ ■ ■ A d+i,T, 

where A m ,T are the barycentric coordinates on T, ej are the canonical Eu¬ 
clidean basis vectors in and oit is a normalizing constant for (p b ,T- The 
function <p b .T is scalar valued and is called a bubble function. The space Qh 
consists of piece-wise linear continuous scalar valued functions. 

Note that if we write v = Vi + v b we have that 


a{u,v) = a(u h vi) + a(u b ,v b ). 


This is so because v b is zero on dT for T e T b and integration by parts shows 
that a(vi,v b ) = 0. We then have the following block form of the discrete 
problem 


( 22 ): 


A 



7 *' 

fi I , where A = 



r<T r^T 



(23) 
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The operators A b , A h G b , Gi and A p correspond to the following bilinear 
forms: 


a(u b , v b ) ->• A b , a(ui, v t ) ->• A h (. KX7p , Vq) ->• A p 
-(divu 6 ,p) = (v b ,Vp) -> -(divv z ,p) ->■ Gg 
v b E V b , u h Vi EVi, p, q E Q h . 


ft is well known that inf-sup condition holds for the MINI element for 
the Stokes problem, and therefore, by Corollary [lj we obtain the following 
inf-sup condition for MINI element discretization of poro-elasticity operator: 
There exists y 0 independent of h, r and K, such that for any (v, q) E V b ,x Q h 
we have 


sup 

(w,s)£V h xQ h 


(A(v,q), (w,s)) 

III M I 


> 7o III (v,q) |||. 


(24) 


As it is well-known (see H), equation (|24| is equivalent to the estimate 


lll(«,p)lll<7o- I |l(/,9)l|. 


(25) 


3.3. Stabilization via elimination of bubbles 

All Pl-Pl stabilized discretizations which we consider here, are derived 
from the MINI element by eliminating locally the bubble functions. For 
details on such stabilizations we refer to the classical paper by Brezzi and 
Pitkaranta [34] (see also 

We now consider the following operator on V) x Q b . 


Ai = 



Gi \ 

-(rAp + St}) ’ 


where S b = G b A b l G b , 


which is obtained after eliminating the equation corresponding to bubble 
functions from (23). This is also an operator of the form given in (16) with 
C = tA p + S b . We have the following theorem: 


Theorem 2. Suppose that the triple ( u b ,ui,p ) solves 



Then the pair ( ui,p ) solves 



(26) 


( 27 ) 
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Moreover, a uniform inf-sup condition such as (24) holds: For any (vi,q) G 
Vi x Qh, 

C Ai(vi,q),(wi,s )) 

SU P - m/.„ „MH -— 7i|(w Q)l\- (28) 

(■ wi,s)eVixQ h 


Proof. Since ( Ub,ui,p ) solves the system (26) we have that 

u b = —A b l G b p 
Aiui + Gip = fi. 

Gju, + G T b u b - t A p p = g — ^ T - ^ T 


G t ui - (G b A b G b + rA p )p = g. 


From this we conclude that ( ui,p ) solves (27). Now, since ( u b ,ui,p ) solves (26), 
from (25), 

IIIK,^,P)||| < 7o" 1 |l(0) fu9)\\i 


and therefore we have 


111(^^)111 < l\(,u b ,ui,p)\\ <'y 0 1 \\{OJi,g)\\ = To II- 


This estimate shows that Ai is a bounded isomorphism, which is equivalent 
to the inf-sup condition ( |28j ). This completes the proof. □ 

Applying Corollary [2] to A/ then shows that any operator C : Of, Q' h , 
spectrally equivalent to rA p + S b will result in a stable discretization of the 
Biot’s model. As we show in the next section (Theorem [3|, the perturbations 
spectrally equivalent to S b are of the form 


(Cp,q) = C T h 2 T 

T&Th 


(Vp- V?), 


where Ct, T G 77 are constants independent of the mesh size h or r. 


3-4. Perturbations, spectrally equivalent to the Schur complement 

In this section we compute the Schur complement (the perturbation or 
the stabilization) given by S b = GjA^Gh. We denote V b ,r = span ip b> r 
and we have that V b = ©Ter h H,T- Let ny be the number of vertices in the 
triangulation, rir be the number of elements, and n b = dnx■ Note that n b 
equals the dimension of V b . With every element T G 77 we associate the 
incidence matrices It G M n vx(d+i) anc j j t ^ ^n b xd ma ppj n g the local degrees 
of freedom on T to the degrees of freedom corresponding to Q and V b . 
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Let us now give a more precise definition of the incidence matrices It 
and Jt for an element T e Th, with vertices (j\ ,..., jd+i), jk £ {1, - - -, ny}, 
and j t j m , for j 7 ^ m. Let {hi, ..., h rf +i}, {ei, ..., e nv }, {/ 1 , .. and 

{Vh ■ ■ ■ iVd} be the canonical Euclidean bases in M d+1 , M nv , M nt and M d , 
respectively. We also denote by (. ki,...,ka ) the degrees of freedom corre¬ 
sponding to the bubble functions associated with T £ Th- We then define 

d -\-1 d 

K-vxH+i, 3 It = E *”**'' 3 * = E (29) 

m =1 ra=l 

Since the sets of degrees of freedom corresponding to the bubble functions 
in different elements do not intersect, we have Jf Jr = Idxd, and, Jf,,Jr = 
0 when T' 7 ^ T. Here Idxd G M. dxd is the identity matrix. Using these 
definitions, we easily find that 

A b = ^ Af l 

T&T h 

G b =Y^ JtG^ t I t t . 

TeT h 

These identities then give, 

S b = GlA~ 1 G b , and hence S b = ^2 I T G bT A~^G b)T lT ■ (30) 

TeT h 

We next state a spectral equivalence result which shows that S b introduces 
a stabilization term of certain order in h for Pl-Pl discretization. Such 
stabilization techniques have been discussed by Verfiirth in |36]J (see also 
§ 8.5.2 and § 8.13.2 in [33]). 

Theorem 3. Let L be the stiffness matrix corresponding to the Laplace op¬ 
erator discretized with piece-wise linear continuous finite elements. Then the 
following spectral equivalence result holds 

S b ~ h 2 L, (31) 


E 

T£T h 


■JTA h rpjrp j 


where the constants hidden in ” are independent of the mesh size. 


Proof. The spectral equivalence is a direct consequence from Lemma [IT] and 
the relations given in (30). D 
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Remark 4. The spectral equivalence in Theorem [3] and the analysis that 
follows justifies the addition of stabilization terms to both the MINI element 
and the stabilized PI-PI discretizations. The results in \Appendix A\ also hold 
for one, two and three spatial dimensions and also give the exact perturbation 
(stabilization) to PI-PI elements that provides inf-sup condition with the 
same constant as the MINI element. 

Related results (in 2D) are found in a paper on Stokes equations by Bank 
and Welfert where it was shown that in 2D the elimination of the bub¬ 
bles in the MINI element gives the Petrov-Galerkin discretization by Hughes, 
Franka and Balestra 138 1 / and Brezzi and Douglas UFA . Here we not only com¬ 
pute the exact Schur complement in any spatial dimension, but we also show 
that the perturbation is spectrally equivalent to a scaling of the discretization 
of the Laplacian with piece-wise linear finite elements. The details are in the 
appendix. 

Such results, however, do not say anything about the monotonicity of 
the corresponding discretization (except in ID, where a further stabilization 
can be introduced in order to obtain a monotone discrete scheme). In fact, 
for the one dimensional case considered in detail in Section the minimum 
amount of stabilization that provides monotone discretization can be calcu¬ 
lated precisely. In general, even for two and three spatial dimensions, adding 
a stabilization term of the form ch 2 L in case when L is a Stieltjes matrix 
improves the monotonicity properties of the resulting discrete problem. This 
is natural to expect because a Stieltjes matrix is monotone. Indeed, the nu¬ 
merical results that we present later also show that adding such stabilizations 
leads to monotone schemes. However, no theoretical results on the mono¬ 
tonicity of the discrete operators for two and three dimensional problems are 
available in the literature and seem to be very hard to establish. 


4. Error estimates for the fully discrete problem 

In this section, we consider the error analysis of the finite element dis¬ 
cretization of the Biot’s model. To simplify the notation and without loss of 
generality in this section we assume that the boundary conditions for both 
the displacement u and the pressure p are homogeneous Dirichlet boundary 
conditions. Then, the weak form of the Biot’s model is as follows: Find 
u(t ) G [iFg(f2)] cZ and p(t) G Hq(UI), such that 

a(u,v) - (div v,p) = (f,v), Vv G [#o(^)] d > ( 32 ) 

-(div d t u, q) - a p (p,q) = 0, Vq E H^(fl), (33) 
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with the initial data u( 0 ) and p( 0 ) given by the solution of the following 
Stokes problem: Find w(0) G [/Jq and p(0) G L 2 (fl), such that, 

a(w(0),u) - (divu,p(0)) = (f(0),v), Vv G [tfo(^)] d , ( 34 ) 

— (divu(0), q) — 0, Vg G L 2 (Q), (35) 


We consider the fully discretized scheme at time t n , n — 1, 2,..as the 
following: Find u\ = u h (t n ) G V h C [ H 1 (^)] <i and p £ = p h (t n ) G ft C 
// 1 (f 2 ), such that, 


a ( u h>v h ) - (divv h ,p%) = (f(t n ),v h ), Vv h eV h , (36) 

-(di vdtujU, q h ) ~ a p (ph, q h ) - eh 2 (yd t pl , \7q h ) = 0, Vq h G (37) 

where d t u% ■— (v% — u %~ 1 )/r and 8 t p^ := (p£ — /4 _1 )/t. Here we try 
to analyze MINI element and stabilized Pl-Pl element in a unified way, 
therefore, the finite element spaces 14 and Qh denote both Stokes pairs. We 
also define the following norm on the finite element spaces: 

/ \ 1/2 

\\{u,p)\\tA-= [\\u\\l + r\\p\\l p +£h 2 \\Vp\\ 2 j . (38) 


We further denote, by || • ||fc and j • 4 the norms and seminorms in the Sobolev 
space and without loss of generality, by || • || the L 2 (f 2 ) norm, i.e. 

|| • || = || ■ ||o- Below we also denote by c a generic constant independent of 
time step, mesh size and other important parameters. 

For the initial data u° h and p° hl we will consider two cases. First case is 
that they are given by the following stabilized Stokes equation: 


a(u° h ,v h ) - (divu h ,p£) = (f(0),v h ) \/v h G V h , (39) 

-(di vu° h: q h ) - £/i 2 (Vp°, Vg h ) = 0 Vq h G Q h . (40) 


Second case is that they do not satisfy (39) and (40) but are defined as 
following, 

di vu° h = 0 and p° h = 0. (41) 


To derive error analysis of the fully discretized scheme (36)-(37), we need 
to define the following elliptic projections Uh and ph for t > 0 as usual, 


a(u h ,v h ) - (divv h ,p h ) = a(u,v h ) - (dwv h ,p), \/v h G V h (42) 
a P (Ph, qh) = a p (p, q h ), \/q h G Q h (43) 
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To estimate the error, following Thomee, [30] we split the discretization 
error as follows. 


u(t) - u h (t ) = (u(t) - u h (t )) - (u h (t) - u h (t )) =: p u - e u , (44) 
p(t) ~ Ph(t ) = (p(t) - p h (t )) - (pr,(t) ~ Ph(t)) —■ P p — e p . (45) 


For t — t n we use the short hand notation p™ = p u (t n ), and similarly e™, p™, 
e™ denote the values of e u , p p and e p at time t = t n , respectively. 

For the error of the elliptic projections, because we use MINI element or 
Pl-Pl element, we have, for all f, 


11 Pu 11 a < ch(|w| 2 + |p|l), 

IIPpIIi < ch\p\ 2 , 11 Pp 11 a p < ch\p\ 2 

\\Pp\\ < ch 2 \p\ 2 . 


(46) 

(47) 

(48) 


We refer to [27j| for details. Since d t p = d t p, we have the estimates above also 
for d t p u and d t p p , where on the right side of the inequalities we have norms 
of d t u and d t p instead of norms of u and p respectively. 

The following lemmas estimate the error between the elliptic projection 
{uh(tn),Ph(t n )} and the numerical solutions {u^,p^}. 


Lemma 5. Letw{ := d t u(tj) — Uh ^ “ h< ' tj ^ andw J p := dtp{tj) — Ph(t ^ P T h(tj 


we have 


n 



(49) 


If the initial data u° h and p® satisfy (39) and (40), we have, 




and if the initial data u° h and p° h are defined by (41) and do not satisfy (39) 
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and (40), we have, 


\e n W < 
|e pl l“ P — 


Vm 


( e !i e p)llT,h + cr 1/2 


1/2 




J 112 
w II a 


d=i 


1/2 


1 / 2 ' 


E eh2 n Vro pii ; 

o'=i 


^£/l 2 ||Va t p(^- 

^'=1 


.( 51 ) 


Moreover, we also have the following estimate in the L 2 -norm, 

n 

ll e pII < c||(e°,eJ)|| Ti/l + cr^ (|K|| a + £ 1 / 2 /i||Vi^|| +£ 1 / 2 /i||Va i p(t i 

3 =1 


(52) 


Proof. Choosing u = Vh G 14 in (32) and q = Qh G £*4 in (33), and subtracting 


both equations from (|36j) and (37), and we have for all Vh G 14 and qh G Q/, 

a«,«h) - ( div Vh, e”) = 0, (53) 

(div <9 t e”, q h ) + a p (e£, g h ) + eh 2 (Vd t ep, X7q h ) 

= (div w™,q h ) +4 2 (Vto",V%) - £h 2 (Vd t p(t n ), Vq h ). (54) 


Choose Vh = d t e™ in (53) and qh = e™ in (|54j) and add these two equations 
together, we have 


eZ,e;)\\U = a(eC,er 1 ) + £/ i 2 (Ve^,Ve^ 1 ) 

+ 7 -(div<, e”) + r£h 2 (VWp, Ve£) 
- T£h 2 (\7d t p(t n ), Ve") 


(55) 


< IKI|a||er 1 || a + £h 2 ||Ve-||||Ve-- 1 | 


+ r|| divtu”||||e”|| + r£h 2 || Vw"|| || Ve"|| + re:h 2 || V%>(£„)|| || Ve" 


Thanks to the inf-sup condition (18), and (53) we have 
(divide”) , 

||ep|| < csup— tt— 7-^- + ci£ 1/2 h\\ Ve"|| 

= c sup + Ci£ 1/2 h|| Vep|| = c||e"|| a + ci£ 1/2 h\\ Ve"||.(56) 

Vh&Vh ||V/i||a 
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Note, for MINI element, we have C\ = 0 and, for Pl-Pl element, c\ > 0. 
Therefore, 


«K)K < 


l<IUI|er 1 || a + ^ 2 ||Ve-||||Ve-- 1 || 

+cr||<||a (HcSIU + Ci£ 1/2 h||Ve-||) 
+r£h 2 ||V<||||VeN|+rf||V^(t n )||||Ve^ 


which implies 


(Ce;)\u< 


K-\e n p - 1 ) 


T, h +cr (\\w™\\ a + £ 1/2 h\\ViUp\\ + £ 1/2 h\\S7d t p(t r 


We sum over all time steps and we have the estimate (49) 
For the error estimate of e™, from ( |53| ), we have, 

a(3 t e™,v h ) - (di vv h ,3 t e*) = 0. 


(57) 


Note that, if the initial data u° h and satisf y (p9| ) and (40), (57) holds for 
n = 1,2,3,.... Otherwise, for initial data (41), (57) only holds for n = 
2,3,... 


Choosing Vh = d t e " in (57) and qh = d t e™ in (54) and adding the two 
equations, and we have 


T- 1 \\p n — P n ~ l II 2 4- ||p n || 2 
' 11 c u 11 a ‘ II 11 a p 


< llp n ll Up” -1 
— 11 c p I1 a P I1 c p 


pll 
71—1 I 


+ T£h 2 1| X7d t e, 
a P + II div<||||e: 

+T£h 2 \\W Wp\\ || V^e”|| + T£h 2 \\W d t p{t n )\\ \\S7d t e, 

KlUKKk + c ll div <ll (IK - < _1 |U + ci^ 1/2 /i||V(e” - e™~ 1 

+reh 2 || VWp || || V9tep|| + T£h 2 \\'Vd t p(t n )\\\\'V d t e, 

1 

\\p"-\r + 

2 


< -||e n || 2 + -He 
— 2 H%> ' o "v 


pi 

n—l || 2 


o n | 

-pi 


= n | 
-pi 


I a 


p 


+ ct\\ div w™\\ 2 + t l \\el 


n- 1 1 | 2 


a + -T£hr\\Vd t e$ 

+\reh 2 \\Vw;\\ 2 + ^T£/i 2 ||V<9 t e£|| 2 + ^T£h 2 || v^p(t n )|| 2 + l -T£h 2 \\Vd t e n p 


where we use the inf-sup condition (56) to estimate ||e 
have 


0 n 

-p 


~n— 11 


Now we 


\K\\lp ^ K iap + c (r|K || 2 + T£h 2 || Vw^ll 2 + T£h 2 \\Vd t p(t n )\\ 2 ) . (58) 

Now we need to consider two different cases due to the initial data. If the 


initial data satisfy (39) and (40), then above inequality (58) holds for n — l 


and by summing up from 1 to n, we can get (50). 
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If the initial data is only defined by (41), (58) does not hold for n — 1 
anymore, we need to estimate 1111 separately. In order to do that, we take 
n — 1 in (55) and then use the inf-sup condition (56) to estimate ||e3,||, 


\ el u\\l + r \\ el p Wa P 


+ ^ 2 ||VeJ || 2 = a(ele° u )+eh^el,Ve° p ) 

+ r(div^,ej) + re/i 2 (VwJ, VeJ) 
- T£/i 2 (V<9 f p(fi), VeJ) 

1„ n„o 1 


< -lie 

— 2 II u\\a 1 2 * I 11 a. 1 2 P 


«„ 0 +olKlla+2^ H Ve pll +2 £/l livin' 

111 1 nO 1 , On — 1 nO 3 


+ cr 2 KH 2 + ^IKII 2 + p^ 2 || Vedl 2 + )Veh 2 || W ,|| 2 


6 

1 \ 211V7 „ 1112 , 3 2_7 2IIV70 MI2 , 1 i,2| 


1112 


+ -eh \\We \\ + -r 2 eh 2 ||V^p(t 1 )|| 2 + -eh 2 ||Ve£|| 


6 


6 


This means 


.bl 2 

-'p 11 a p 


< ^rll(e°, e° p )\\ T:h + c (t||u3 || 2 + re/r 2 || VwJ || 2 + re/r 2 || V<9 t p(ti)|| 2 ) . 


Now we summing up (58) from 2 to n and use above estimate of ||ep|| 2 p , we 
can get (51). 


Finally, the estimate (52) follows directly from (49) and (56). 
Next lemma give the estimations of w J u and wL 


□ 


Lemma 6. Letu(t ) andp(t ) be the solution of (32) and (33), w J u = d t u(tj ) — 
an dp u (t) = u(t)—u h (t). Assume d tt u(t) G L 1 ((0, T], [Ffg(r2)] d )n 
L 2 ((0,T],[i?£(fi)] d ) and d tt p(t) G ^((0, T], H^Q)) D L 2 ((0, T], we 

have, 


E 

3 = 1 
n 

E 

3 = 1 


KIU < c 


0 


Zn 1 /*^n 

|| <?««!! idt + - J ||9ip u ||idt ) , 


\w: 


j ||2 ^ 



c r 


||a tt ti|| 2 df + - / ||<9 t p u || 2 df 


r 


(59) 

( 60 ) 


Moreover, let w J p = d t p(tj ) 


Ph{tj)-Ph(tj- 1) 


and p p = p{t ) — Ph{t). we have 


E» v bii 


< c 


i=i 


1 Z*^ 71 


^l|V ^|| 2 <c(r 
t=i 



’itPllidf + - 

r 


ttp nidi + - 


ll^iPpllldt , 


ll<%> P ||idt 


(61) 

(62) 
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Proof. We consider 


l( J — j Q tU (l-^ — | ( u ifj) u {tj-l) Uhitj) Uh(tj- 1 ) 


Note that 


then we have 


r 


= : <l+^,2- 


K,i = - 

T 


1 f tj 


[s - t j ^ 1 )d tt u(s)ds, 


'U- 


■j—1 


1 f tj 

K,2 = - dtpu(s)ds , 




KIU < IK^IU + IK^ 

-ii r 


1 [ tj 

(s-£j_i)<9 tt u(s)ds|| a +-|| / 9 t p u (s)ds|| a 


< c 


'tj -1 
ftj 

'tj-1 


T 


tj—1 


1 r tj 


||9 tt «||ids + - / ||||ids , 


r 




j -1 


then (59) follows directly. Moreover, we have 

1/2 


All? < 


1 / 2 ' 


rl/2 / ll^lli ds +r 1/2 / ||0 t p u ||ids 


'U- 


j-1 


'ti~ 


3-1 


1 


< c[t ||^ t n||jds+ - / ||%>Jids , 


'tj-i 


T 


'tj -1 


then (60) follows directly. Estimates (61) and (62) can be obtained similarly, 
which completes the proof □ 


Assuming extra regularities of the exact solutions u(t) and p(t) as usual 
for convergence analysis of the finite element method, we have the following 
theorem about the error estimates for the error (u — Uh)(t n ) and (p — Ph)(t n ). 

We assume that u and p have all the regularity required by the proof of 
the theorem below, which more precisely means that, for q = 1, 2, oo and 
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s = 1, 2 we have: 


u(t) € L” ((0,7’], [ff 0 W) nl“ ((0,7], ]/7 2 (fi)] d ), 

d,u(t) € L‘ ((0,7], (J7 2 (!J)]‘ i ) , d u u(t) e L* ((0,7], [ff 0 1 (O)]‘ i ) , 

P(<) 6 7” ((o, 7], 7J‘(S1)) n L” ((0,7], ff 2 (Ji)), 

9.p(«) 6 £"(( 0 ,7], 7/^(!Ti)) n i“(( 0 ,7], J7 2 (n)), d u p(t) 6 L* (( 0 ,7], 7/ 2 (n)) 


Theorem 7. Lei w(t) and p(t) be the solution of (32) and ( [33]) 7 u^ and p 1 ^ 
be the solution of @ and (|37|). For displacement u{t), we have 


( u(t n ) u hiP(^n) Ph ) Ik,/! 


< II ( e °> e ?) \\r,h+C<T 



+h 


||0 tt u||idt + / e 1/2 h\d tt p\idt 
Jo 

rtn 

u(t n )\ 2 + \p(t n )\ 1 +(r 1/2 + £ 1/2 h)\p(t n )\ 2 + / (\d t u\ 2 + \d t ph) dt 


e 1/2 h\d t p\ 2 dt 


+ t n max £ 1//2 h\\\7d t p(tj)\\ 


l<j<n 


(63) 


For pore pressure p(t), if the initial data u® and p ° satisfy (39) and (40), we 
have, 


\Wn)-pl\\a P 


< ||e p || a p +c{T 


in \ 1/2 

IMI?diJ 


in \ 1 / 2 ' 


+h 


|p(^n) 12 + 


in \ 1/2 / /*in \ l/ 2 

(|^w| 2 + |9tp|i) 2 diJ + i J £h 2 \d t p\ldt) 


+Vtn max 

l<j<n 





(64) 
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If the initial data u° h and p® are defined by (41), we have 

\\p(t n ) ~ Ph\\a p 
1 


< 




K4) 


r.h 


+ c < r 


\ V2 

\\d tt u\\ldt\ + 


+h 


\pifn) 12 + 


(\d t u\ 2 + \dtp\x) 2 dt 


1/2 


+ 


rtn \ 1 / 2 ' 

£h 2 \\d a p\\{dt 


\ V 2 ' 


+V^ max £ l l 2 h\\Vd t p(tj)\\ > . (65) 

i <j<n J 

Moreover, for pore pressure, we also have the following error estimate in 
L 2 -norm, 

I \p(tn) ~Ph\\ 

( T fitn rtn 

||9 tt M||id/+ / £ 1/2 h\d tt p\idt 


< c ll { e °u, e p) II T,h + c\r 


+h 2 \p(t n )\ 2 + h 


I "IL 

(|0 t u| 2 + |^p|i)dt+ / £ 1/2 h\d t p\ 2 dt 


Uo 


+t n max £ 1/2 h\\S7d t p{tj)\\ 

l<j<n 


( 66 ) 


Proof. The estimate ([63]) follows directly from (|44j), (45j), ( |46[ ), ( 47]) , 

(59), (61), and triangle inequality. Note that we used (46) and (47) not only 
for u, p, but also their counterparts for d t p u and d t p p . 


Similarly, ( |64[ ) follows from ( |45| ) ( |5d| ), ( [60] ) ( [62] ) ( |46| ) ( |47[ ), and their 
versions for the time derivatives of the error and the triangle inequality. 


Next, for the second set of initial conditions, (65) follows from (45), (51), 
(60), (62), (46), (47) (applied also for time derivatives of the error), and the 
triangle inequality. 

Finally, (66) follows from (45), (48), (52), (59), (61) and the triangle 
inequality. □ 

Remark 8. All the error estimates in Theorem^ consist of two parts. One 
part is the error for t > 0 which, in all cases, gives optimal convergence 
order. The other part is the error in the approximation of the initial data, 
i.e., 


°, e p)llnh an d \\e° p \\a p . From the triangle inequality, we have 

[ e u,e° P )\\ T ,h < c[\\(p 0 u ,p° p )\\ Tjh + || (u(0) -u° h ,p(0) -pi) || Tifc ] , 


-'pWUp 


< HA P + lb( 0 )-P°Ja P , 
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where p° u and p p are the errors due to the elliptic projection and (it(0) — w°) 
and (p(0) — p°) are the errors due to the choice of initial conditions, either 
satisfying Stokes equation (39) and (40) or the simpler given in (41). 

If the initial data satisfies the stabilized Stokes equation (39) and (40), 
the initial errors strongly depend on the regularity of the initial data. A cru¬ 
cial role is played by the assumptions on the regularity of the pore pressure 
p(0). If we assume p(0) G Hq(Q), then the standard error estimates for the 
elliptic projection and stabilized Stokes equation show that the initial data 
errors are appropriately bounded, and, hence, we have optimal order of con¬ 
vergence for the discrete scheme. Therefore, the overall convergence rate of 
the stabilized MINI element is optimal. However, if we assume that p(0) is 
merely in L 2 (fA), then we cannot expect that the errors in the initial data are 
of optimal order, and, therefore, the overall convergence rate of the stabilized 
MINI element is not optimal as well. 

If we just use the simple practical choice (41), we cannot expect thatu° h p° h 
approximate w(0) and p(0) in general. Therefore, regardless of the regularity 
assumption of the initial data, the overall convergence rate of the stabilized 
MINI element will not be as desired. However, in some cases, even when 
the initial errors are large, they decay with respect to time (see As a 

consequence, the discretization error when using stabilized mini element is 
still optimal for sufficiently large time (long time). 


5. Numerical Experiments 

In this section, we present several numerical experiments in order to il¬ 
lustrate the performance of the proposed stabilized methods. We will choose 
well-known benchmark problems in order to deal with different aspects as 
variable permeability, different boundary conditions, the accuracy of the ap¬ 
proximations, etc. 


5.1. Layered porous medium with variable permeability 

In the first experiment we want to illustrate non-monotone pressure be¬ 
havior when we have a low permeability in a sub-domain. We consider 
a test proposed in m which models a porous material on which a low- 
permeable layer (K = 10~ 8 ) is placed between two layers with unit perme¬ 
ability (K = 1), as shown in Figure 5.1 The boundary of the square domain 
is split in two disjoint subsets Tj and T 2 on which we assume the following 
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Figure 5.1: Domain representing a square of layered porous material with 
different permeability. 


boundary conditions: on the top, which is free to drain, a uniform load is 
applied, that is, 

p = 0, a ■ n = g, with g = (0, —1)*, on r 1; (67) 

whereas at the sides and bottom that are rigid the boundary is considered 
to be impermeable , that is, 

Vp ■ n — 0, u = 0, on r 2 . (68) 


Zero initial conditions are considered for both variables, and the time step is 
chosen as r = 1. Notice that this test can be reduced to a one-dimensional 
problem. Therefore, in the following simulations we will show the numerical 
solutions corresponding to one vertical line in the domain as displayed in 
Figure |5.1 


First we approximate using linear finite elements for displacements and 
pressure. If no stabilization term is added to the discrete formulation, the 
approximation for the pressure field that is obtained by using 32 elements 
on the grid is shown in Figure 5.2 (a). We observe that strong spurious 


oscillations appear in the part corresponding to the low-permeable layer. 
However, if the stabilized scheme is used for the simulation with the same 
number of nodes, the oscillations are completely eliminated and the method 


gives rise to the monotone solution for the pressure, as we see in Figure 5.2 
(b). 
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(a) (b) 

Figure 5.2: Numerical solution by Pl-Pl for the pressure to the two-material 
problem (a) without stabilization term and (b) with stabilization term. 


Next, we use approximation by MINI element with the same number of 
elements. Similarly to the previous case, when no stabilization parameter is 
included in the formulation, the oscillatory behaviour of the pressure approx¬ 
imation is evident, as shown in Figure [5~3| (a). Notice that the oscillations are 
much smaller than in the case of Pl-Pl elements, but are still not eliminated 
by using this Stokes stable pair of spaces. Again, a perturbation stabilizes 
the method and we obtain oscillation-free approximation for the pressure 
held (see Figure 5.3 (b)). 


5.2. Mandel’s problem 

Mandel’s problem (see [32]) is an important benchmark problem because 
the analytical solution in two dimensions on a finite domain is known. It is 
an excellent model that can be used to verify the accuracy of a discretization. 
Mandel’s problem models an infinitely long poroelastic slab sandwiched at 
the top and the bottom by two rigid frictionless and impermeable plates. 
The material is assumed incompressible and saturated with a single-phase 
incompressible fluid. Both plates are loaded by a constant vertical force 
as shown in Figure |5.4, where a 2a x 2b wide cross-section is displayed. 


The force of magnitude 2 F per unit length is suddenly applied at t — 0, 
generating an instantaneous overpressure by the Skempton effect [33], which 
will dissipate near the side edges as time progresses due to the drainage 
effect, since the side surfaces {x = ±a) are drained and traction-free. In this 
problem, it turns out that the horizontal displacement u is independent of 
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Figure 5.3: Numerical solution by P2-P1 for the pressure to the two-material 
problem (a) without stabilization term and (b) with stabilization term. 


the vertical direction y, whereas the vertical displacement v is independent 
of the horizontal coordinate x. The analytical solution for the pore pressure 
can be found in j44j and is given as follows 


p(x,y,t) = 2*p 0 Y^ 


sm a r . 


n= 1 


Oir, 


sm a„ cos a r 


a n x 


cos 


cosa r 


exp 


-a^ct 


(69) 

where po = ^B( 1 + u u )F, being B the Skempton’s coefficient that for our 
problem is B — 1 and u u = the undrained Poisson’s ratio, c is 

the consolidation coefficient given by c = K (A + 2/i), and a n are the positive 
roots of the nonlinear equation 


tan otr, — 


1 - V 


CYr, 


- V 


As can be observed in (69), also the pressure is independent of the vertical 

shows that the normalized pressure is 


direction. In fact, Coussy (see 
the solution of the following equation 


dp d 2 p 
dt dx 2 


71=1 


a 2 sin a n cos a n . 9 ^ 

——— -— exp (-alt) 

a n — sm a n cos a n 


(70) 


Note that the right-hand side is constant in space and it can become large 
at the beginning of the process. 
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Figure 5.4: 2D physical and computational domains for Mandel’s problem. 


For the finite element solution, the symmetry of the problem allows us 
to choose only a quarter of the physical domain as a computational domain, 
as shown in Figure |5.4 Moreover, the rigid plate condition is enforced by 
adding constrained equations such that vertical displacements on the top are 
equal to an unknown constant value. The triangulation of the computational 
domain is obtained from a uniform rectangular grid n x x n y by splitting each 
element in half. The dimension of the porous slab is specified by a = b = 1, 
and the material properties are given by K = 10 -6 , E = 10 4 , v = 0, and 
therefore v u = 0.5. The Lame coefficients are computed in terms of the 
Young modulus and the Poisson ratio as follows, 


A = 


Ev 


E 


(1 -2z/)(l + !/)’ 


h = 


2(l + i/)' 


Finally, the applied force has a magnitude of F = 1M Pam. 

The first test with Mandel’s problem will illustrate the need of stabilizing 
the Pl-Pl discretization, as well as the MINI element discretization, in order 
to remove the spurious oscillations in the pressure field. We choose a final 
time T = 10“ 4 for the computations with only one time-step, and a spatial 
grid with n x = n y = 32. Since the pressure unknown is independent of the 
vertical coordinate, we will present the results on a representative horizontal 
line. In Figures 5.5 (a) and 5.5 (b), we show the numerical solution for the 


pressure (plotted in circular symbols) obtained by using Pl-Pl finite element 
methods without and with stabilization, respectively. The numerical solution 
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Figure 5.5: Numerical solution by Pl-Pl of the pressure for Mandel’s problem 
(a) without stabilization term and (b) with stabilization term. 


is plotted against the analytical solution that is displayed by a dashed line. 


The same comparison is shown in Figures 5.6 (a) and 5.6 (b) for the MINI 


element scheme. For the latter, the inf-sup condition is satisfied, but we ob¬ 
serve that nonphysical oscillations appear in the pressure field, albeit smaller 
than in the Pl-Pl case. By adding in both methods stabilization terms, 
oscillation-free solutions are obtained, as seen in Figures [5.5[ b) and 5.6 [b). 
Next, we analyze the behavior of the pressure in different times. For this 
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Figure 5.6: Numerical solution by MINI element of the pressure for Mandel’s 
problem (a) without stabilization term and (b) with stabilization term. 


purpose, in Figure 5.7 the solution of the pressure obtained by stabilized 
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Pl-Pl finite elements on a grid with n x = n y = 32, together with the corre¬ 
sponding analytical solution are shown in different times. We can observe a 
good agreement between both solutions for all the cases. A very interesting 
behavior of the solution of Mandel’s problem is that it can achieve values 
greater than one at some time instants. In the literature, this is known as 
the Mandel-Cryer effect and usually is associated to a lack of monotonicity. 
However, it is clear that this phenomenon is due to the source term that ap¬ 
pears in equation (70), and is fully in agreement with the maximum principle 
for the heat equation. 



Figure 5.7: Comparison of numerical and analytical solutions of the pore 
pressure for Mandel’s problem at various times. 


Finally, we investigate the convergence properties of the proposed sta¬ 


bilized schemes by comparing the analytical solution, given in (69), with 


the numerical solution obtained on progressively refined computational grids 
with n x = n y ranging from 10 to 80 and with time-steps (r = T/n t ) from 0.5 
to 0.0625. In Table 5.1, for each mesh and a final time of T — 1, we display 


the error for the pressure in the norm 

\\p(tn)-Ph \\ 2 = Mtn) - Phil 2 + KT\\S 7 {p{t n ) ~ PhW 


Frorn Table 5T we observe first order convergence, according to the error 
estimate obtained in Theorem [7} A very interesting insight rising from these 
results is that similar errors for both finite element methods are obtained. 
This is due to the fact that very similar stabilization parameters have to 
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n x x 7i y x n t 

10 x 10 x 2 

20 x 20 x 4 

40 x 40 x 8 

80 x 80 x 16 

Pl-Pl 

0.0163 

0.0110 

0.0058 

0.0029 

MINI 

0.0162 

0.0110 

0.0058 

0.0030 


Tabic 5.1: Energy norm of the error for the pore pressure by using 
stabilized Pl-Pl and MINI element for different spatial-temporal grids. 


be added to both methods to avoid the nonphysical oscillations, since the 
addition of the bubble plays a positive role but with a very small contribution. 
This point could be a reason to support the use of the stabilized Pl-Pl 
scheme against the MINI element that also has to be stabilized. 

5.3. Barry & Mercer's problem 

Another well-known benchmark problem on a finite two-dimensional do¬ 
main is Barry & Mercer’s model, see m It models the behavior of a rectan¬ 
gular uniform porous material with a pulsating point source, drained on all 
sides, and on which zero tangential displacements are assumed on the whole 
boundary. The point-source corresponds to a sine wave on the rectangular 



0 317 a 

V = u = —= 0 
dy 


Figure 5.8: Computational domain and boundary conditions for the Barry 
and Mercer’s source problem. 


domain [0,a] x [0,6] and is given as follows 

f{t) = 2f35 (X0M) sm(f3t), 


(71) 
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where f3 = 


(A + 2fi)K 
a b 


and 5( X0i y 0 ) is the Dirac delta at the point (xo,yo). 
In Figure [578] the computational domain together with the boundary condi¬ 
tions are depicted. The boundary conditions do not correspond to a realistic 
physical situation, but they admit an analytical solution making this model 
a suitable test for numerical codes. Here we use this model to assess the 
monotone behavior of the approximations of the pressure. 



oo o 0.5 1 


(b) T = 3tt/2 

Figure 5.9: Numerical solution for the pressure by Pl-Pl and deformation of 
the grid after applying the pulsating pressure point source, for two different 
values of T. 


We consider the rectangular domain (0,1) x (0,1), and the following values 
of the material parameters are considered E = 10 5 , v = 0.1 and K = 10~ 2 . 
The source is positioned at the point (1/4,1/4) and a right triangular grid 
with n x = riy = 64 is used for the simulations. The solution for the pressure 
produced by the stabilized Pl-Pl scheme is plotted in Figure |5.9| for two 
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different “normalized times” t = fit of values t = 7t/2 and t = 3ir/2. Also we 
display the deformation of the considered triangular grid, according to the 
results obtained for the displacements. We can observe that depending on the 
sign of the source term (positive for t = 7t/2 and negative for t = 3n/2) the 
resultant displacements cause an expansion or a contraction of the medium. 

The analytical solution of this problem is given by an infinite series, and 
can be found in m- It has been observed that solutions displayed in Fig¬ 
ure 


5.9 resemble the exact solution very precisely. 


Fluid pressure oscillations for the Barry and Mercer’s problem can be 
demonstrated by considering the standard schemes given by a Pl-Pl or MINI 
element discretizations. In order to see this characteristic non-physical os¬ 
cillatory behavior, a small permeability and/or a short time intervals are 
considered. Therefore, in the previous test, we have changed the value of K 
to 1CT 6 and T to 10~ 4 . For these parameters, in Figure 5.10 we show the 
numerical solutions obtained for the pressure held, by using Pl-Pl scheme 
(on the top) and the MINI element (on the bottom). We can observe that if 
no stabilization term is added to any of the discrete schemes (left pictures), 
then non-physical oscillations appear in the surroundings of the source-point. 
However, by adding the proposed artificial stabilizations, we can see (right 
pictures) that these oscillations are completely eliminated. 


6. Conclusions 

In this paper we have analyzed the convergence and the monotonicity 
properties of low order discretizations of the Biot’s consolidation model in 
poromechanics. While the convergence results are complete in some sense, 
there are still several open theoretical questions regarding the monotonic¬ 
ity of the resulting discretizations. Clearly, our numerical results show that 
choosing the stabilization parameters correctly lead to oscillation-free solu¬ 
tions, but justifying this rigorously is difficult and a topic of ongoing research. 
We have to say though that as a rule of thumb, one can choose stabilizations 
that are optimal in ID, and, the resulting approximations in higher spatial 
dimensions will be oscillation-free. 
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Figure 5.10: Numerical solution for the pressure field by Pl-Pl (top) and 
MINI element (bottom), without and with stabilization term, at a final time 
of 1(T 4 and a permeability of K = 10~ 6 . 
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Appendix A. Local elimination of bubbles 

In this appendix we compute the contribution of bubble stabilization 
in the MINI element. We show that G^A^Gb is spectrally equivalent to 
the stiffness matrix corresponding to the discretization of the Laplace with 
continuous piece-wise linear finite elements. 

To begin, we fix T e Th and we prove several simple identities. When 
the dependence on T neds to be emphasized we indicate this by indexing the 
corresponding quantities with T, but most of the time, this is not needed 
and we set 


Afc = A k,T, k — 1,..., (d + 1), 

a = ay, and ip = ip hT = ccAi... A^+i. 

Here A k,r(x) are the standard barycentric coordinates on T and is a 
constant chosen so that (pb,T has a value 1 at the barycenter of T. To integrate 
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polynomials over a d-dimensional simplex T we use the well known formula 
for integrating powers of the barycentric coordinates (see my 



\T\ 


ffi! • ■■Pd+M\ 

{(5\ + ... + j3 d+ 1 + d)\ 


(A.l) 


Further, we introduce the matrix A E M dx ( d+1 ) whose columns are the ap¬ 
propriately scaled gradients of A*,, k = 1,..., (d + 1) i.e. 


A= vlTf(VA 1 ,...,VA d+ 1 ) = VW\ 


/VAi • ei, 

\VAi • e d , 


VArf + i • e\ 

VAd + i • e d . 


We note that A T A equals the local stiffness matrix for the Laplace equation 
on T, namely 

{Lr)jk = (A t A )jk = j VAfc • VAj. 


Note that we have 


d -\-1 d -\-1 

Vy? = aJ^XfcVAfc, Xk = A j, k = 1,..., (d + 1). 

k =l j=i;j¥=k 

With this notation in hand, we now prove two auxiliary identities. 
Lemma 9. For V(p we have 

(i) j(\7<p\7ip T ) = a 2 r] d AA T , 

(ii) J T \V(p \ 2 = a 2 ri d tr (L T ). 

2 d ~ l d\ 


Here r) d = 


(3d)! ' 


Proof. To prove (i) we observe that = — j=A\, X — (An 

vm 

Since A is a constant matrix (independent of x) we have that 




a 


A ( l XX 1 ) AA 
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The formula given in (A.l) gives that 


T 


, T XX ') i rJr XiXk = ,,dm u, 


2, j = k, 


Hence, J T XX T = Vd\T\(I + 11 T ), where 1 = (1,..., 1) T . As A fc = 1, 


we have that, Ylt=i VA k = 0, or, equivalently, A1 = 0. These identities show 
that 

j = a 2 r] d A(I + 11 T )A T = a 2 ri d AA T , 

and the proof of (i) is complete. 

To show that (ii) holds we observe that f T |V(/?| 2 = J T tr(Vy?Vy9 T ), and 
we can use (i) to compute that 


d-\-1 


/ |V</?| = / tr(V^Vv? ) = tr I / Vy?V<^ ) = a 2 r] d tr(L T ). 
Jt Jt \Jt 

In the last step we used that tr(AA T ) = tr(A T A). 


□' 


Using this lemma we now calculate the local stiffness matrices for Ab t T 
and Gb,T- 

Lemma 10. For Ab,T and Gb.r we have 


r b,T 


(i) A btT = a 2 r] d (iitr(L T )I + (A + fi)AA T ). 

(n> Gt - T ~ (mTI)! 

Proof. To show the identity for A^t recall that 

( \r)jk = a(((pe k ),((pej)) 

= 2/r / e(((pe k )) : e((<pej)) + A / div(^e fc ) div(^e i ). 


A straightforward calculation shows that 

1 


= ~(A/^e k + e k (V^) ), 
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and hence 


5,- 


e((m)) : £ ((^ e j)) = ir / |V<A + ^( / 




We also have 


div(<y3e fc ) div(^ej) = / S/tpV(p 


jk 

Finally, using Lemma [9] for A,t we get 

At = a 2 77 d (Aitr(L T )/ + (A + /j)AA t ). 

To show (ii), we have, for k — 1,..., (d + 1) and j — 1,..., d, 

A ik f 


(fib^r)jk / i}P^j ' ^Afc) (V A/j • 6j) IP 


vTl 




(A.2) 


Computing / T ^r concludes the proof of (ii). 

From this, for the local Schur complement S^t we get 

S b , T = Gj T A-^,T = c d |T|A T (/itr(L T )/+(A + /i)AA T )' 1 A 
= aA T (/ + /3AA t )~ 1 A 


□ 


where 

c d \T\ X + fi _ d\(3d)\ 

' “/itr (LrY P ~fitr(L T y ° d ~ 2 d -\(2d + l)!) 2 ' 

We apply the Sherman-Morrison-Woodbury to obtain that 

(/ + M T A) _1 = / - PA t (I + f3AA T )~ 1 A 


This then shows that 

aA T (I + /3AA t )- 1 A = [I - (/ + /3A)- 1 ] . 

Observing that 

/ - (/ + At) -1 = /-((/ + At) - At)(/ + At) -1 = At(/ + At) - 1 , 
we obtain that 

Sfj^T = <?L t {I + At) (A.3) 

We next show that S^t behaves as a scaling of the local stiffness matrix 
corresponding to the Laplace operator. 
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Lemma 11. We have the following spectral equivalence result 

Sb,T ~ h^Lx, 

with constants independent of the mesh size. 

Proof. Let /+ > /x 2 > ... > Hd > /++1 = 0 be the eigenvalues of the scaled 
matrix L T = tr( ] L , r) Lt, and ifi,... ,ifd+ 1 be the corresponding eigenvectors. 
Note that because of the scaling, wejrave that p,k can be bounded indepen¬ 
dently of the mesh size hr- We set f3 = /3tr (Lx) = ~~ and we obtain the 

following representation of Lx'- 

d c? 

l t = y^ ( j + plx )- 1 = {i+n T r= 

j =i j= 1i +ppj 


Obviously, similar relation holds for S^x because the eigenvectors of Lx (Lx) 


and Sb,T are the same (this is easily seen from (A.3)). We then have that for 
any x G M. d the following inequalities hold 


(Sb,xx,x)p — 


c d \T\ 


E 




)l 7=i 1 + fil-L 




Hence, 


Cd\T\ 


/i(l + ( 3 / j.i 


■( L t x,x)p < (S b ,TX,x)p < 


c*\T\ 




(L T x,x)p. 


We write everything in terms of Lx, and from the obvious relations tr(L^) fa 
h^j7 2 , \T\ fa hx we conclude the proof of the lemma. □ 

Remark 12. As is easily seen, for d = 1 we have that both bounds coincide, 
and in fact, we have that 


Sb,T — 


Mi + P) 


-Lt — 


hf 


12(2p +A)' 


(A.4) 


where we have used that /i\ — 1 and \T\ = hx in Id. 
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